Timing analysis techniques at large core distances for multi-TeV gamma ray 

astronomy 



V. Stamatescu' , G. P. Rowell, J. Denman, R. W. Clay, B. R. Dawson, A. G. K. Smith, 
T. Sudholz, G. J. Thornton, N. Wild 

School of Chemistry and Physics, The University of Adelaide, Adelaide 5005, Australia 
t now at Institut de Fi'sica d'Altes Energies, Edifici Cn., E-08193, Bellaterra (Barcelona), Spain 



Abstract 

We present an analysis technique that uses the timing information of Cherenkov images from extensive air showers 
(EAS). Our emphasis is on distant, or large core distance y-ray induced showers at multi-TeV energies. Specifically, 
combining pixel timing information with an improved direction reconstruction algorithm, leads to improvements in 
angular and core resolution as large as ~ 40% and ~ 30%, respectively, when compared with the same algorithm 
without the use of timing. Above 10 TeV, this results in an angular resolution approaching 0.05°, together with a 
core resolution better than ~ 15 m. The off-axis post-cut y-ray acceptance is energy dependent and its full width 
at half maximum ranges from 4° to 8°. For shower directions that are up to ~ 6° off-axis, the angular resolution 
achieved by using timing information is comparable, around 100 TeV, to the on-axis angular resolution. The telescope 
specifications and layout we describe here are geared towards energies above 10 TeV. However, the methods can in 
principle be applied to other energies, given suitable telescope parameters. The 5-telescope cell investigated in this 
study could initially pave the way for a larger array of sparsely spaced telescopes in an efibrt to push the collection 
area to > 10 km^. These results highlight the potential of a 'sparse array' approach in efl'ectively opening up the 
energy range above 10 TeV. 

Keywords: gamma-ray astronomy, cosmic rays, imaging atmospheric Cherenkov telescopes, instrumentation 



1. Introduction 

Based on current experimental evidence ob- 
tained from instruments such as H.E.S.S., MAGIC, 
CANGAROO-III, VERITAS and MILAGRO, it is 
generally accepted that particles can be accelerated in 
shell-type supernova remnants to energies beyond a 
few xlO'"^ eV. The H.E.S.S. instrument has observed 
a number of_y-ray sources with strong fluxes above 
10 TeV m |2> H, with some sources showing no 
evidence of a spectral cut-off. Furthermore, most of 
the new sources found on the galactic plane have hard 
differential energy spectra with F < 2.5 f?]. Future 
observations of ~ 100 TeV y-ray emission from these 
acceleration regions will conclusively test current 
acceleration models and determine if SNR and/or 
other types of sources are responsible for cosmic ray 
acceleration up the knee \^^\M- 
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An important advantage of the multi-TeV regime 
(>10 TeV) is that it becomes easier to distinguish be- 
tween the hadronic (y-rays from cosmic -ray interac- 
tions) and leptonic (y-rays from the inverse-Compton 
or IC process on soft photon fields) emission chan- 
nels. Electrons suffer strong energy losses due to syn- 
chrotron cooling if strong magnetic fields are present in 
their acceleration region. At the same time, the Klein- 
Nishina effect reduces the IC cross section. These ef- 
fects strongly suppress the IC process so, unless the y- 
ray emission is associated with continuous or sporadic 
sources of multi-TeV electrons (i.e. pulsar wind nebu- 
lae, active galactic nuclei), observed spectra from steady 
and extended galactic sources that reach into the tens of 
TeV may favour a hadronic scenario. 

The current generation of Imaging Atmospheric 
Cherenkov Telescope (lACT) arrays are limited bytheir 
collecting area in the multi-TeV energy range P, [sl- 
Good spectral measurements at energies approaching 
100 TeV and beyond, assuming reasonable observation 
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times (~ 50 h), may require a dedicated instrument with 
collecting areas of around 10 km^ l"?]. Based on expe- 
rience with H.E.S.S, a wide field of view is also de- 
sirable, both in terms of its survey potential and for 
the purpose of selecting off'-source regions for cosmic- 
ray background determination in imaging large scale 
sources. Meeting these design goals while maintain- 
ing a good angular resolution is important to the astro- 
physics potential of any such future instrument. This 
is a key performance parameter that directly impacts on 
the point-source sensitivity and on the ability to probe 
the morphology of extended sources. 

In this paper we investigate the use of shower tim- 
ing information to improve stereoscopic shower recon- 
struction. This follows from our earlier work \2j and is 
motivated by our 'sparse array' design fl, with which 
showers are often viewed at large core distances (typ- 
ically > 200 m) by a cell of five telescopes, and, also 
by the first simulations of a multi-Te V cell by lUt] . The 
term 'core distance' refers to the distance between the 
telescope and the shower axis, calculated in the shower 
plane. In an alternative approach, Colin et al. [9'] con- 
sidered the performance of a dense lACT array over the 
energy range of 1 to 100 TeV. Large-scale telescope ar- 
rays of varying densities, aimed at covering an energy 
range up to ~100 TeV are also under consideration for 
the planned CTA observatory [10.1 . 

The fluorescence technique Ml ill pioneered air shower 
reconstruction methods of detectors such as Fly's Eye 
and HiRes, in which angular and temporal informa- 
tion from a single shower image was combined to re- 
construct the shower geometry. For the atmospheric 
Cherenkov imaging technique, the longitudinal and 
transverse temporal structures of EAS were first inves- 
tigated using the HEGRA I ACT system 1 12]. By con- 
sidering the image pixels times, it was found that a so- 
called time gradient exists along the image major axis, 
and that it is correlated with the core distance. At small 
core distances (< 100 m) this correlation is not present 
for hadronic showers, and this difl'erence, as suggested 
in lUl], is used by MAGIC to improve back- 

ground rejection. This approach was also investigated 
using simulations of a wide field of view (10°) camera 
IitIi . using a single telescope as well as a 2-telescope 
system. In that paper it was shown that the time gradi- 
ent could be used to improve background rejection by 
a factor of 2 - 3 for the case of a single telescope. Al- 
though we do not attempted to link angular resolution 
directly with background rejection in this study, we note 
that incorporating our timing analysis in a more compre- 
hensive multi-parameter analysis, such as that used by 
iITtIi . merits further investigation. 



The time gradient is also used by VERITAS to de- 
fine an optimum signal integration window [13], which 
improves the Cherenkov signal to skynoise ratio of the 
recorded signals. Timing information was used for 
time-based image cleaning by the CANGAROO col- 
laboration (e.g. 1 18j |19J) in order reduce the effect of 
night sky background. It is also used by MAGIC as a 
way of lowering the standard image cleaning thresholds 
HE. 



2. Simulations of a 5-telescope cell for multi-Te V en- 
ergies 

Our simulations of EAS use the CORSIKA v6.204 
1 2111 and SYBYLL li221 packages, coupled with a tele- 
scope response simulation based on simJelarray ||20[ 
2311 . Gamma-ray showers are generated over the energy 
range 1 -500 TeV using a flat dN/dE oc E" energy spec- 
trum to enable good statistics at high energies, and with 
a zenith angle of 30°. The chosen observation altitude 
of 220 m a.s.l. is typical of many Australian sites. 

The simulated cefl layout, shown in Figure [1] is such 
that four telescopes are situated in the corners of a 
square of side 500 m, while the fifth telescope is at its 
centre. Each simulated shower is used 20 times by scat- 
tering the core location over a 1 km radius (in the shower 
plane) with respect to the centre of the lACT cell. Sim- 
ulations of the cell performance in which showers were 
re-scattered 10 times gave consistent results, in terms 
of angular resolution as well as the mean and r.m.s. of 
width and length, when compared to our default (20 
times re-scatter) method. We note the use of 10 times re- 
scattering of CORSIKA showers was successfully em- 
ployed in the analysis of real H.E.S.S. data ll24ll . The 
collection area achieved with this cell layout approaches 
1 km^ above 10 TeV, and exceeds Ikm^ above 100 TeV 
I2I. 

Each simulated telescope has //1. 5 optics with a 
modest-sized 6 m diameter mirror of elliptic profile 
lEill . of area 23.8 m^ (comprising 84 x 60 cm spherical 
mirror facets). The circular camera has an 8.2° diam- 
eter field of view (FoV) made from 804 square pixels 
of side length 0.24°. Previous simulation studies 1127 1 
found that a pixel diameter of ~ 0.25° allowed for a 
good quality determination of the shower image sec- 
ond order moments (see discussion below), while opti- 
mum image orientation parameters were achieved using 
~ 0.15° diameter pixels. In the present study, the choice 
of pixel size is motivated by the size of the optical point 
spread function of each telescope. Ray-tracing results 
li2&l indicate that the elliptic dish shape provides a well- 
centred off-axis blur spot and < 0.25° focussing out to 
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Figure 1 : Simulated cell layout: four telescopes are at 
the comers of a cell of side length L - 500 m with the 
fifth telescope at the centre. 



~ 4° off-axis. Compared to a traditional Davies-Cotton 
design, this dish shape gives a small improvement of 
~ 10% in the 4° - 5° off-axis angle range. 

A MODTRAN 1 29] tropical atmospheric profile and 
maritime haze model are used to treat the wavelength- 
dependent attenuation of Cherenkov light by the atmo- 
sphere. In addition to the atmospheric transmission, 
simJelarray treats the detection of photons by mod- 
elling the mirror reflectivity, the photomultiplier (PMT) 
collection efficiency and the photo-cathode quantum ef- 
ficiency, as discussed in 1I20I1 . The resulting photoelec- 
tron (pe) signal is then input into a detailed electron- 
ics simulation that models the data acquisition and tele- 
scope trigger 

Single pe pulses (after the PMT, pre-amplifier and 
shaping amplifier) are added in time over a 100 ns buffer 
such that the median Cherenkov pe arrival time is ~ 40 
ns after the start of the memory buffer The result- 
ing signal in each data acquisition channel, which is 
AC-coupled and includes fluctuations from night sky 
background (NSB) and electronic noise, is continuously 
sampled at 1 GHz and digitized using a 12 bit ADC 
(analog to digital converter). An analog memory buffer 
would allow for this data acquisition scheme to be im- 
plemented in hardware. A low gain channel is used if, 
at any given time, the digitized signal of the high gain 
channel exceeds the dynamic range of the ADC. 

The pixel signal values and their corresponding pulse 
arrival times are obtained using a sliding integration 
window. The highest digitized value in the time buffer 
is found, the digitized signal in a 21 ns region around 
the peak is integrated and the number of pe in the chan- 
nel is determined. This integration time is sufficiently 



long to capture a single PMT pulse for the range of core 
distances relevant to this study. The time of any given 
pixel (with respect to the start of the time buffer) is given 
by the time bin that corresponds to the digitized signal 
peak. 

Single pe fast pulses are added in the time domain to 
test if a channel triggers its discriminator These pro- 
duce output signals that are subsequently summed to 
determine whether a camera has triggered. A trigger 
signal for each telescope results if the signal in 3 near- 
neighbour (grouping of 9) pixels exceeds 6 pe, and the 
centre pixel in the grouping is compulsory. This con- 
dition is based on simulations ll30ll that indicated a low 
accidental single telescope trigger rate (~10"^ Hz) due 
to NSB, and the need to accept as many events as pos- 
sible in order to obtain a large effective area. A soft- 
ware stereo trigger is required for an event to be recon- 
structed. 



A standard two-level image cleaning method Bill is 
used to reduce the effect of NSB on shower images. 
Pixels with signal over the (higher) image (or picture) 
threshold are included in the image, while those over the 
(lower) boundary threshold must have a next neighbour 
above the image threshold to be included. A previous 
investigation of the effect of image cleaning on the cell 
performance 1I25I1 suggested a choice of 8 pe image and 
4 pe boundary thresholds. 

Shower images are parametrized using the Hillas for- 
malism II32II . which yields the Hillas parameters: a ma- 
jor axis, an image centre of gravity (cog), and shape pa- 
rameters width and length, which are second order mo- 
ments in the direction perpendicular and parallel to the 
major axis, respectively. The parameters prescribe an 
image as an ellipse. The image size parameter is the 
sum of the pe in the pixels that pass cleaning. 

Given that a shower image is a two dimensional pro- 
jection of the shower, the major axis is an approxi- 
mate projection of the shower axis in the camera plane. 
Therefore, if the image is well parametrized, the ma- 
jor axis will pass close to the true shower direction in 
the camera plane. The pointing direction of the major 
axis is not prescribed by the Hillas formalism. This axis 
'pointing degeneracy' may be broken if another major 
axis from a second telescope is available, since, for well 
parametrized images, the axes will intersect in the vicin- 
ity of the shower direction in a common camera plane 

Shower reconstruction is performed for events in 
which at least 2 telescopes trigger and pass the follow- 
ing quality cuts: an image size cut of > 60 pe and a dis2 
cut of < 4°. dis2 is the maximum modulus value the of 
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the vector sum: 

dis2 — maxQcog + axis length]) 



y (deg) A 



(1) 



where cog is the angular distance from the camera cen- 
tre to the cog, and axis is the unit vector of the image 
major axis. The cut on dis2 mitigates the camera edge 
effects on shower images, and further optimization of 
this this cut is a topic of a future study. 

3. Shower reconstruction 

Using our standard event reconstruction, the shower 
direction is determined by a weighted mean of the ge- 
ometric intersection points of the major axes in a com- 
mon camera frame. This scheme, which was developed 
for the HEGRA lACT system, is denoted Algorithm 1 
in ll34ll . The weights used in calculating the mean of 
intersection points are given by w: 



(sizCa + sizet) sm(s) 



(2) 



where a and b denote two images used to obtain the 
intersection point, and s is the so-called stereo angle be- 
tween the two major axes. This gives more weighting 
to brighter image pairs, which will, as a result of bet- 
ter Hillas parameterization, have more accurate major 
axes. This choice also gives less weight to major axis 
pairs that are nearly parallel, which results in poor di- 
rection reconstruction fsT]. The weights w do not influ- 
ence the reconstruction if only one pair of major axes is 
available (i.e. = 2). Once the reconstructed shower 
direction is found, the core position is determined using 
a similar method, by intersecting major axes in the re- 
constructed shower plane, starting from the positions of 
the telescopes . 

We also employ an improved method of reconstruc- 
tion known as Algorithm 3 |34], which was also devel- 
oped for HEGRA and which has been used in H.E.S.S. 
analysis (e.g. BSll ). Algorithm 3 is illustrated in Fig- 
ure |2] for the case of only one pair of major axes. For 
each shower image, the expected uncertainties in the 
major axis orientation and the cog position are used to 
define an error elUpse. These uncertainties are obtained 
from lookup tables with dependencies on image size and 
width/length. The lookup tables are generated from sim- 
ulations, using a flat dN/dE oc E" spectrum, by consid- 
ering the positional dififerences between the parameter- 
ized major axis of each image and the true shower axis 
projected into the camera plane [28]. This gives an es- 
timate of the true source position together with errors 
associated with the estimate. The error ellipses are then 
combined using a weighted mean (e.g. see 1361 281'). 
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Figure 2: The case of only two available images in an 
event reconstruction. The Algorithm 1 direction recon- 
struction is given by the intersection of the two major 
axes. Two error ellipses are combined using Algorithm 
3 to give an improved estimate (X = (x, y)) of the true 
source position. The predicted distance dp for each im- 
age is used to prescribe the angular displacement of 
each error ellipse from its respective cog, along each 
major axis. The true distance, d,, is the angular distance 
between the true source position and the image cog. 



where the weights are chosen to minimize the error el- 
lipse associated with the Algorithm 3 result (X = (x, y) 
shown in Figure |2]i. The errors used to prescribe each 
each error ellipse are then re-determined by using this 
improved source estimate and the weighted mean is re- 
calculated. The Algorithm 3 source position usually 
converges to a stable point after two or three iterations. 

The quantity dp, which is illustrated in Figure |2] is 
used to prescribe the position of each error ellipse on 
its corresponding major axis. It is the expected angular 
distance between the image cog and the true source po- 
sition, and is traditionally predicted using combinations 
of image parameters, such as length or width/length cou- 
pled with image size. The accuracy and precision with 
which dp is predicted is a contributing factor to the per- 
formance of Algorithm 3. 

4. The Cherenkov image time gradient 

The temporal and angular structure of an air shower 
may be sampled by using the pixel times and their an- 
gular camera plane coordinates, respectively. Figure |3] 
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Figure 3: Time-profiles of five 30 TeV y-ray showers 
simulated from zenith and imaged at several core dis- 
tances, which are indicated on each panel. Individual 
time-profiles are denoted by different colors. The pixel 
time t is with respect to the weighted (using pixel ampli- 
tude) mean of image pixel times, whilst dist is the pro- 
jected distance along the major axis of each pixel with 
respect to the cog. 



shows the time-profiles of five example 30 TeV y-ray 
showers. For each shower image, it shows the pixel 
time, t, versus dist, the angular distance of each pixel 
from the cog, projected along the major axis. Here the f- 
axis origin corresponds to the pixel amplitude-weighted 
mean of recorded pixel times. Figure |3] clearly shows 
that the average slope of individual time-profiles tends 
to increase with a larger distance between the telescope 
and the shower core. 

The EAS time development detected by lACTs was 
investigated by Hess et al. 1 12]. It was found that when 
telescopes are located within the Cherenkov shoulder 
(~ 150 m iIstIi ). the detected photons that are emitted 
by ultra-relativistic secondary particles low in the atmo- 
sphere tend to arrive early in time compared to those 
emitted near the top of the shower In our study the fo- 
cus is on y-ray showers detected at large core distances. 
For showers detected beyond ~ 150 m from a telescope, 
photons emitted near the top of the shower tend to ar- 
rive before those emitted near the bottom. We have in- 
vestigated this effect for showers up to 10° off-axis, a 
Hmit which is imposed by our field of view. Showers 
viewed at large core distances are also more extended 
in time, which can be seen in Figure [3] The time gradi- 
ent parameter is determined for each Cherenkov image 
as the average slope of its time-profile. In our analy- 
sis this is calculated using a weighted linear fit, where 
each weight is the relative signal strength (i.e. number 



of pe relative to the image size) of an image picture pixel 
(with > 8 pe). 

We use a three dimensional toy model of y-ray 
shower longitudinal time development to examine the 
correlation between time gradient and core distance, 
and the dependence of time gradient on the height of 
maximum Cherenkov light emission. It can also be used 
to test the effect of off-axis shower geometries, which 
indicates that the time gradient remains approximately 
constant for a shower axis that is inclined with respect to 
the telescope axis by angles of < 10° (l^. The model is 
best applied to showers with core distance greater than 
about 150 m and provides a check on the simulations in 
the large core regime (> 200 m from a telescope). 

The toy model handles the case of a vertically pointed 
telescope, illustrated in Figureffl Using an approach 
similar to that of Cabot et al. [38], we calculate the time 
delay dt at the telescope, for a given core distance x: 

dt(x,z) = - (L-z/cos(y)) (3) 
c 

where dt is the difference between the arrival time at the 
telescope of light emitted on the shower axis at height z 
and the time at which the extrapolated primary y-ray 
reaches the ground level, if it were propagating with 
speed c. The optical path L in Eq. [3] is obtained by 
integrating the altitude-dependent refractive index rjih) 
along the path of the the Cherenkov light. A simple de- 
pendence of the refractive index on altitude is assumed: 
T](h) = 1 pjpe-''^''", with T]o = 2.76 ■ lO^^ and ho = 8.0 
km (see 02811 for details). Using angles y and 6 as de- 
fined in Figure |4] (a), this leads to an expression for the 
time delay given by Eq. lA.ll in [Appendix A| 

Figure |4] (b) shows that the light reaching the tele- 
scope is imaged into the camera plane to a point with 
angular coordinates (px,Py), given by: 

Px(x,z) - cosx py(x, z) - sinx (4) 

where the angles Q. and x determined by Eq. IA.2I 
and Eq. IA.3I respectively in |Appendix A| 

The results of the toy model are compared in Figure 
|5] with Monte Carlo simulated time profiles obtained 
for an on-axis 30 TeV y-ray shower Differences be- 
tween the model and the simulations become evident at 
smaller core distances (< 100 m) where the transverse 
development of the shower, which is not modelled, be- 
comes important. For the toy model, the origins of the 
t and dist axes correspond to the height of maximum 
light emission, Zmax- 

Zmax 180 X 
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Figure 4: The geometry of time delay and imaging toy model: (a) A shower axis inclined by angle y to the telescope 
axis and with azimuth angle 6 lands at core distance x from the telescope. Cherenkov light (shown in red) emitted at 
height z on the shower axis, reaches the telescope inclined with respect to the vertically pointed telescope axis by O 
and with azimuth angle x- (b) Light whose direction is prescribed by the angles D. and is imaged into the camera 
plane to a point (p^,pv)- 



where c/, is the angular distance between the true source 
position in the camera and the image cog, which is illus- 
trated in Figure |2] In choosing z - Zmax as the reference 
point, the dist-axis origin in Figure |5] corresponds, by 
definition, to the cog. For the f-axis, dtizmax, x) is only 
approximately equal to the weighted mean time of im- 
age pixels (w.rt. the time of the extrapolated primary 
reaching the ground). This assumption does not hold as 
well for showers with small core distances, as seen from 
the the negative f-axis-ofFset of the curve with respect 
to the simulated points in Figure |5] (a). However, the 
choice of axis origin does not alter the modelled time- 
profile shape, which is what determines the toy model 
time gradient: 



Mt _ dtiz + dz) - dtiz) 
AQ ~ Q(z -1- dz) - Q.{z) 



(6) 



We choose dz - 0.1 km, so that the quantity Adt/AQ. 
may be considered a local time gradient, which corre- 
sponds to a part of the shower around a given value of 
z. 

Hess et al. fl?] found that the time gradient depends 
strongly on the distance between the telescope and the 
shower core. This is also apparent from Figure which 
shows the correlation between these quantities. The 
'reflected' points in Figure |6] with positive time gra- 
dient at core distances > 200 m, are caused by events 



with poorly reconstructed Algorithm 1 source positions, 
which are used to break the pointing degeneracy and re- 
sult in the major axis pointing away from the true source 
position. The data in Figure|6]are color-coded according 
to height of maximum Cherenkov light emission, which 
is calculated for each image as: 



h, 



t/,7r/180. 



cos(zen) - r sin(zen) 



(7) 



where dr is the angular distance in the camera between 
the cog and Algorithm 1 reconstructed source position, 
and where r is the Algorithm 1 reconstructed core dis- 
tance (in the reconstructed shower plane). 

In order to compare Ac/f/AQ with with our simula- 
tion results, we adjust the toy model to handle telescope 
pointing at non-zero zenith angles. For simplicity, only 
on-axis shower geometries are considered and this is il- 
lustrated in Figure|7] Time and angular development are 
modelled in terms of height z and core distance x by Eq. 



A.4l and Eq. IA.5I respectively in Appendix A 



The efltect of hi on the time gradient is modelled by 
choosing z and z+dz in Eq. |6]that match the correspond- 
ing range of hi used to split up the simulated data. These 
ranges are given in Figure |6](b), together with a distri- 
bution of hi for our y-ray data set. Figure |8] shows that 
the toy model predictions (Adt/AQ.) and simulations are 
consistent for core distances beyond ~200 m. It also 
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Figure 5: Simulated pixel time t (black points), with respect to the weighted mean time of image pixels, versus (dist), 
its distance along the major axis from the cog, for an on-axis 30 TeV y-ray shower from zenith with core distance 
from the telescope of 80 m (left), 200m (middle), 400m (right). Solid lines are the results of the toy model with angles 
7 and 6 set to 0, and z ranges from 1.5 km to 15.0 km (left), 2.5 km to 15.0 km (middle), 5.0 km to 15.0 km (right). 
The section of each line with the most positive values of dist corresponds to the top of the shower axis. 
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Figure 6: The scatter plot of time gradient versus true 
core distance, using the simulated y-rays with primary 
energy between 1 TeV and 500 TeV. The major axis de- 
generacy is broken using Algorithm 1 (see text for de- 
tails). The colors correspond to the indicated ranges of 
reconstructed light maximum, hi. 
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Figure 7: Modified time delay and imaging toy model 
incorporating a zenith angle but simplified to handle 
only on-axis geometry. The x-axis is the line connecting 
the telescope to the shower core, with the telescope po- 
sition taken as the origin. The true core distance is taken 
as X cos(zen), the distance between the telescope and 
the shower core projected onto the true shower plane. 
Light (in red) emitted at height z from the shower axis 
reaches the telescope inclined with respect to the tele- 
scope axis by D., which is measured as an angular dis- 
tance in the camera. 
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Figure 8: The time gradient calculated from simulated 
data (points) and toy model (solid lines) are plotted 
against true core distance projected in the shower plane. 
The matching zto z + dz range in the model and hi range 
used to select data are color coded according to Figure 
|6l In the case of simulated data, the major axis pointing 
degeneracy, which affects the sign of the time gradient, 
is broken using the true source position in the camera 
(this only applies to the data presented in this plot). 



shows that in the large core distance regime the depen- 
dence of time gradient on h/ is weak. The toy model 
time gradient is not applicable for small values of x, 
where shower images become tend to be rounded and 
have poorly defined major axes. 

Given that the angular distance d, will increase with 
core distance, a strong correlation also exists between 
di and the time gradient. This is evident in Figure|9]and 
suggests that the time gradient may be used to predict 
the angular distance J,. From Figure |9]it is also apparent 
that the accuracy of this prediction can be improved if 
the dependence of c/, on the parameter/!/ is taken into ac- 
count. This is because a shower whose maximum emis- 
sion is deeper in the atmosphere will result in a larger 
value of df. Conversely a shower with the same core 
distance that has its maximum emission at a higher alti- 
tude will result in a smaller dt- 

A lookup table with dependencies on time gradient 
and 1 //?/ is filled with average values of d,. This table is 
then used in the analysis to obtain the predicted angular 
distance dp. As detailed in Section [3] this angular dis- 
tance prescribes the position of a predicted source posi- 
tion along each major axis. The predicted uncertainties, 
which are obtained using the error lookup tables, define 
an error ellipse around this predicted source position. In 
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Figure 9: The scatter plot of angular distance d, versus 
time gradient, using the simulated y-rays with primary 
energy between 1 TeV and 500 TeV. The colors corre- 
spond to the ranges of reconstructed light maximum, hi, 
given in Figure |6] 



this way, shower timing information is used in the Al- 
gorithm 3 direction reconstruction. 

It worth noting that the optimal way in which to use 
pixel timing information may be as part of a multi- 
parameter gamma-ray maximum likelihood reconstruc- 



tion (e.g. rt39L 14011 ). rather than the piece-wise approach 



employed here. However we leave this for future work. 

5. Cell performance using Algorithm 3 

5.7. Predicting the distance to the source 

In this study we compare the performance of Algo- 
rithm 3 using timing information with results obtained 
using traditional predictors of the distance d,. The two 
lookup tables used for this comparison have dependen- 
cies on length and logio(image size) and on width/length 
and logioiimage size), respectively. Figure [TO]shows the 
behaviour of dt in terms of the dependencies of these 
lookup tables. 

The length of an image increases with core dis- 
tance, which makes it a predictor of li,. A dependency 
on logioiimage size) is included in this lookup table 
because the image length increases with image size, 
as shown in Figure [10] (a). The image aspect ratio, 
width/length, decreases with increasing core distance as 
shown in Figure [10] (b). Both predictors are affected by 
camera edge effects for the distance range of d, > 3°, 
which results in otherwise elongated images becoming 
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Figure 10: The true distance dr (deg.) plotted versus the predictors (a) length and (b) width/length. The colors denote 
ranges of logio(image size) given in (c). 



more rounded. This effect leads to smaller length and 
larger width/length values. 

Both length and width/length suffer from a similar 
problem in predicting dt, which can be seen in Figures 
[TOl(a) and[TOl(b): beyond dt values of ~ 1.2°, each pre- 
dictor is insensitive to changes in dt- Moreover, splitting 
the length, width/length, or length /log\o(miage size) ac- 
cording to hi does not solve this problem due to the 
inter-dependent way in which hi and core distance influ- 
ence the length, width, as opposed to the time gradient 

Beyond di of ~ 2.75° the predictors are well behaved, 
but it is necessary to distinguish between this regime, 
and that of d, < 1.2°. The width/length table is therefore 
split into three separate lookup tables that correspond to 
three ranges of d,: < 1.25°, between 1.25° and 2.75°, 
and > 2.75°. The length lookup table is split in two 
ranges of d, < 2° and dt > 2°. The choice of which 
lookup table to use for a given image is made based on 
the value d,-, which is initially obtained from the Algo- 
rithm 1 reconstruction. 

Images with a small c/, usually correspond to show- 
ers with small core distances, whose time profiles tend 
to suffer from fluctuations [12). We therefore make 
an additional refinement when using the time gradient 
and hi lookup table: for images with d,- < 0.85, the 
width/length distance predictor is used instead. 

The performance of the three sets of lookup tables 
dependencies is illustrated by Figure [TT] The breaks in 
the scatter plots correspond to the aforementioned split 
ranges of dt. The combination of time gradient and hi 
produces the strongest correlation (by eye) between the 



true distance, dt, and the predicted distance, dp. 

5.2. Angular resolution 

We evaluate the angular resolution of the telescope 
cell using r68, the angular radius containing 68% 
of events. We apply image shape cuts on mean- 
scaled parameters 114 ill , which are used for cosmic 
ray background suppression. The images that corre- 
spond to events passing the cuts are, on average, bet- 
ter parametrized, and this leads to an improved angular 
resolution. The cut parameters are mean-scaled-width 
(MSW), mean-scaled-length (MSL) and mean-scaled- 
Npix (MSNPix) MSNPix is a new parameter, cal- 
culated by scaling the number of pixels in each image 
by the expected number of pixels (obtained from sim- 
ulations, for a given image size and core distance) and 
then averaging over the telescopes in the reconstruction. 
The cut on MSNPix provides additional y-hadron sep- 
aration power (after MSW and MSL cuts) for energies 
above 50 TeV, however it does not significantly affect 
the angular resolution. 

The post-cut (MSW< 1.05, MSL< 1.2, MSNPix< 
1.1) angular resolution obtained with the 'time gradi- 
ent dependent' Algorithm 3, which uses the time gradi- 
ent and Ijhi lookup table, is compared in Figure [T2] to 
that obtained obtained by using length or width/length 
together with \og\Q{image size). For on-axis showers, 
the next-best results are achieved with the length and 
log[()(image size) lookup table, and we term this the 
'length-dependent' Algorithm 3. The improvement in 
angular resolution given by the time gradient dependent 
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Figure 11: Predicted distance dp (deg.) plotted against true distance dt (deg.). dp is predicted using a distance lookup 
table based on (a) time gradient and l/hi, and width/length for dr < 0.85°, (b) length and logioiimage size), and (c) 
width/length and logio(image size). 



Algorithm 3 over the length-dependent Algorithm 3 re- 
construction ranges from ~10% to ~40% in the energy 
range of ~5 to ~200 TeV. The error bars are determined 
on the basis of Poisson fluctuations in the number of 
events (68% of events) used to calculate the radii for 
each energy bin. 

Figure [12] also shows, for comparison, the result of 
the purely geometric Algorithm 1 reconstruction. The 
upturn in the Algorithm 1 angular resolution seen above 
~30 TeV is due to a worsening in the reconstruc- 
tion of low telescope multiplicity (2 or 3 telescope) 
events whose core distance from the array centre is large 
(greater than ~300 m). This subset of events suffers 
from smaller stereo angles (nearly parallel major axes), 
thereby increasing the error in major axis intersection 
points. The Algorithm 3 direction reconstruction miti- 
gates this problem. For instance, the level of improve- 
ment provided by /e«^f/i-dependent Algorithm 3 over 
the purely geometric reconstruction, is ~ 10% below 
100 TeV, and ~ 55% above that energy. Given its better 
distance prediction, the time gradient dependent Algo- 
rithm 3 improves things further In this case, the level of 
improvement over Algorithm 1 is ~ 15% below 10 TeV, 
~ 40% in the 10 to 100 TeV energy range, and ~ 75% 
above 100 TeV. 

The simulation results indicate that by using timing 
information together with the Algorithm 3 direction re- 
construction, an r68 angular resolution of less than 0.1° 
is obtained above ~ 10 TeV, and this value approaches 
0.05° at the highest energies. This confirms the finding 
of an earlier study QTJ, namely that a large collecting 
area can be achieved at multi-Te V energies using a small 



number of telescopes, while maintaining good angular 
resolution. Using our time-dependent stereoscopic anal- 
ysis, we have demonstrated this for a 'sparse array' of 
lACTs [6]. In the present study we show this to be case 
for off-axis performance (see the discussion at the end 
of this section). 

5.3. Core and energy reconstruction 

The improved direction reconstruction provided by 
the time gradient dependent Algorithm 3 can be used to 
improve the shower core reconstruction. A new recon- 
structed core position is obtained by intersecting axes 
defined using the time gradient dependent Algorithm 3 
reconstructed source position and the cog of each image 
passing quality cuts. These newly defined axes are inter- 
sected in the reconstructed shower plane, starting from 
their telescope positions. This time gradient dependent 
'Algorithm 3 core' resolution is compared in Figure [13] 
to that obtained using the /engf/i-dependent Algorithm 
3 reconstruction. With respect to the length-dependent 
Algorithm 3 core resoltuion, an improvement of ~ 20% 
is achieved for energies above a few TeV through the 
introduction of timing information. The purely geomet- 
ric Algorithm 1 core reconstruction is also shown for 
comparison. 

Our shower energy reconstruction is performed us- 
ing a y-ray filled lookup table with dependencies on im- 
age size and the reconstructed core distance to each tele- 
scope. The events used to fill the lookup table are first 
re-weighted in order to smooth out effects of the large 
discontinuity in our Monte Carlo statistics at 10 TeV 
and 100 TeV. An energy estimate is obtained for each 
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Figure 12: Comparison of the angular resolution r68 (radius for 68 % event containment) after all shape cuts (MSW< 
1.05, MSL< 1.2, MSNPix< 1.1), obtained using the time gradient dependent Algorithm 3 implementation (red 
triangles), which uses predicted distance dp based on time gradient and l//i/ and two different versions of Algorithm 3 
(blue squares) that use alternative methods for obtaining dp: (a) dp is obtained using length coupled with logio{image 
size), and (b) dp is obtained using width/length coupled with logioiimage size). The Algorithm 1 results (black stars) 
are shown in both panels for comparison. 
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Figure 13: Core Resolution r68 (radius for 68 % event 
containment) as a function of true energy after shape 
cuts: MSW< 1.05, MSL< 1.2, MSNPix< 1.1. The 
results shown as black stars use the Algorithm 1 core 
reconstruction. Those shown as red triangles use the 
time gradient dependent Algorithm 3 core reconstruc- 
tion, while those shown as blue squares use the length- 
dependent Algorithm 3 core reconstruction. 



telescope, and the estimates from telescopes used in the 
shower reconstruction are combined to give a global en- 
ergy estimate, is^co- This simple approach does not ac- 
count for the effect of the height of shower maximum 
on the reconstructed energy, that is demonstrated to im- 
prove the accuracy of the reconstructed energy |f42|. 

The energy resolution and the energy reconstruction 
bias, obtained using the two Algorithm 3 core recon- 
struction methods, are compared in Figure [14] We de- 
fine 'energy resolution' , r68£, as the radius that contains 
68% of events in the (Erea, - E,rue) IE, rue distribution at 
each true energy bin, where E,n,e is the true energy. The 
error bars are based on the Poisson error in the number 
of events within the 68% containment radius. Both Al- 
gorithm 3 implementations lead to a ~ 35% better en- 
ergy resolution above 100 TeV, compared to the purely 
geometric core reconstruction. The gain in core reso- 
lution given by time gradient dependent Algorithm 3 
over the /engf/i-dependent Algorithm 3 does not appear 
to improve the overall energy resolution. However, fur- 
ther work on the energy reconstruction algorithm may 
be beneficial in exploiting the timing information. 

The term 'energy reconstruction bias' in Figure [14] 
refers to the mean of the {Ereco - E, rue) IE, me distribu- 
tion. The error bars given in that figure are calculated 
as r6?>El V^, where is the number of events in the 
underlying distribution. Thus r68£/ ^In is analogous 
to the 'standard error in the mean' in each E,rue bin. 
The prominent positive energy reconstruction bias seen 



11 



at a few TeV for all reconstruction methods is caused 



by threshold effects (e.g. see 114211 ). which are the re- 
sult of upward fluctuations in pixel values from showers 
with energies just below the threshold. This bias can in 
principle be corrected, however we do not apply such 
a correction in this study, since the effect on our y-ray 
results will be negligible above a threshold of 3 - 4 TeV. 

5.4. Off-axis angular resolution 

The large (4.1° radius) FoV of the simulated tele- 
scopes allows the detection of showers with a large 
range of off-axis angles (out to ~ 7°). The energy- 
dependent off-axis camera acceptance is shown in Fig- 
ure[T5] Acceptances are calculated for each off-axis an- 
gle as the ratio of off-axis and on-axis post-cut (MSW< 

1.05, MSL< 1.2, MSNPix< 1.1) events, using the Al- 
gorithm 3 core reconstruction in the determination of 
mean-scaled cut parameters. Off-axis geometries are 
obtained by decreasing the telescope axis altitude an- 
gle in increments of 1°, starting with 60°, which is the 
on-axis value. 

The off-axis acceptance in Figure [15] has a full width 
at half maximum (FWHM) that ranges from ~ 4° to 
~ 8°, depending on the energy. Given the broadness 
of the off-axis acceptance, it is important to achieve 
a good off-axis angular resolution. Increasing the tilt 
of the telescopes allows high energy showers with pro- 
gressively higher core distances to pass event selection 
cuts. It was therefore necessary to increase the COR- 
SIKA simulated 'throw radius', within which showers 
are scattered, from 1 km to 2 km for primary energies 
above 100 TeV, so as not to artificially limit the number 
of detected events. Distant off-axis showers result in a 
higher proportion of small and/or edge-affected images 
that are not entirely removed by the quality cuts. Such 
images may be poorly parametrized and this, together 
with lower average telescope multiplicities, on average, 
tends to worsen our direction reconstruction. 

In Figure[T6] we compare the off-axis angular resolu- 
tion achieved using the time gradient dependent Algo- 
rithm 3 with that obtained using the length-dependent 
Algorithm 3. We also show the the results of the Algo- 
rithm 1 shower reconstruction. The energy threshold of 
the cell rises with increasing off-axis angle, which leads 
to reduced event statistics at energies close to threshold 
(causing larger error bars), rather than any large change 
in the inherent angular resolution of the cell. The gen- 
eral degradation of angular resolution r68 above 100 
TeV, which is present in all reconstruction algorithms, 
is due to the previously mentioned camera edge effects. 

Figure[T6]shows that above ~10 TeV, the time gradi- 
ent dependent Algorithm 3 gives a larger improvement 



over the length-dependent Algorithm 3 than in the on- 
axis case. While the time gradient is insensitive to in- 
creasing off-axis angles, the behaviour of dt as a func- 
tion of length becomes more degenerate (see ll28ll ). As 
a result, the distance resolution obtained with length 
and logioiimage size) degrades rapidly beyond ~ 2° 
off-axis, making this reconstruction method less robust. 
The results of the next-best Algorithm 3 in Figure[T6]are 
obtained using distance lookup tables that are divided 
into two ranges of c/,: < 4° and > 4°. This was chosen 
'by eye' from the behaviuor of d, as a function of length 
at 3° off-axis [28J. Adjusting this value (between ~ 3° 
and ~ 5°) for each off-aixs angle may improve the re- 
sults, although any differences are expected to be within 
10%. 



6. Conclusion 

We have presented a shower reconstruction technique 
for imaging atmospheric Cherenkov telescopes (lACT), 
which uses pixel timing information. Our aim was to 
improve the reconstruction of extensive air showers im- 
aged at large distances (> 200 m) from the telescopes of 
a sparse array. The technique was investigated using a 
simulated representative array of five modest-sized (6 m 
diameter), wide field of view (8.2° diameter) telescopes, 
which would achieve a large y-ray collecting area for 
the multi-TeV energy regime. 

The longitudinal air shower time development gives 
rise to an image time gradient, which strongly depends 
on the distance between a telescope and the position of 
the shower core on the ground. This time gradient was 
used to predict the angular distance between the image 
cog and the source position in the camera. This pre- 
dicted distance was used in an improved reconstruction 
algorithm ll34ll (also known as Algorithm 3). 

The Algorithm 3 performance was compared to al- 
ternative versions that do not use pixel timing informa- 
tion using y-ray showers simulated in the energy range 
of 1-500 TeV. A noticeable improvement in the angu- 
lar resolution (r68) of the cell was obtained by using 
the time gradient in the reconstruction. The level of im- 
provement over the 'next-best version' of Algorithm 3 
varied between ~10% and ~40% for on-axis showers, 
depending on the energy. 

The angular resolution obtained with Algorithm 3 
was also compared to that obtained with a purely ge- 
ometric reconstruction method, known as Algorithm 1 
1 3411 . The improvement with respect to Algorithm 1 
ranged between ~ 10% and ~ 55% (depedning on the 
energy) before the introduction of timing information. 
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Figure 14: Energy reconstruction bias (left panel) and energy resolution (right panel) of the (Erm, - E,rue)IEtrue 
distribution (see text for details), as a function of the true energy, after all shape cuts: MSW< 1.05, MSL< 1.2, 
MSNPix< 1.1. Ereco is the reconstructed energy and Emte is the true energy. Black stars denote the Algorithm 1 event 
reconstruction, blue squares denote the /e«gf/i-dependent Algorithm 3 event reconstruction and red triangles denote 
the time gradient dependent Algorithm 3 event reconstruction. For clarity, in the left panel, the red triangles are shifted 
higher up in energy by 5%. 
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Figure 16: Angular resolution r68 (radius for 68 % event containment) as a function of true energy and off-axis angle. 
The r68 values are given after all shape cuts (MSW< 1.05, MSL< 1.2, MSNPix< 1.1). Algorithm 1, time gradient 
dependent Algorithm 3 and /e«gf/j-dependent Algorithm 3 results are given as black stars, red triangles and blue 
squares, respectively. 
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and between ~ 15% and ~ 75% following its introduc- 
tion. 

For energies in excess of 100 TeV, an r68 angular res- 
olution approaching 0.05° was achieved with the time 
gradient dependent Algorithm 3, while the core resolu- 
tion was better than 10 m. Using the same method, an 
r68 angular resolution better than 0.1° was obtained out 
to ~ 6° off-axis for energies around 100 TeV. This ap- 
proach may improve the survey capability and extended 
source analysis of future instruments such as TenTen |0] 
and CTA flff). 

While the level of improvement in array performance 
does depend on the energy range of interest, on the tele- 
scope parameters, and the analysis method employed, 
we believe our results highlight the potential of using 
timing information in stereoscopic reconstruction in a 
sparse (~ 500 m spacing) telescope array. Thus our 
study demonstrates the potential of large core distance 
Cherenkov imaging for multi-TeV astronomy, with an- 
gular and energy resolution performance similar to that 
of H.E.S.S. in its respective energy regime. 
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Appendix A. 

The following mathematical expressions relate to toy 
models for y-ray Cherenkov image time profiles, which 
are discussed in Section |4] 

dt{x, z) = — ( ^jx^ +7} sec^ 7 + 2xz tan 7 cos (5 



cos 7 



ho 



yjx^ + z^sec^y + 2xz tany cos 6riQ{\ - e ''''»)) (A.l) 



Q.{x,z) = arctan 



X(^,z) = arctan 



^Jx^ + z^ tan- 7 + 2xz tan 7 cos 6 



^z tan 7 sin S 
x + z tan 7 cos 6 



(A.2) 
(A.3) 



dt{x,z) = -( V?^^^<^^l^n(>en>)2 - 



cos(zen) 



— y/z^ + (x + ztan{zen)yr]o{l - e"'"'°)) (A.4) 
z 



0.(x,z) = arctan(- 



X cos(zen) 
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